{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Thermal behavior of a wall"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Modeling and Simulation in Python*\n",
    "\n",
    "Copyright 2021 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-nc-sa/4.0/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# install Pint if necessary\n",
    "\n",
    "try:\n",
    "    import pint\n",
    "except ImportError:\n",
    "    !pip install pint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# download modsim.py if necessary\n",
    "\n",
    "from os.path import basename, exists\n",
    "\n",
    "def download(url):\n",
    "    filename = basename(url)\n",
    "    if not exists(filename):\n",
    "        from urllib.request import urlretrieve\n",
    "        local, _ = urlretrieve(url, filename)\n",
    "        print('Downloaded ' + local)\n",
    "    \n",
    "download('https://raw.githubusercontent.com/AllenDowney/' +\n",
    "         'ModSimPy/master/modsim.py')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import functions from modsim\n",
    "\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This case study is based on Gori, Marincioni, Biddulph, Elwell, \"Inferring the thermal resistance and effective thermal mass distribution of a wall from in situ measurements to characterise heat transfer at both the interior and exterior surfaces\", *Energy and Buildings*, Volume 135, 15 January 2017, Pages 398-409, [which I downloaded here](https://www.sciencedirect.com/science/article/pii/S0378778816313056).\n",
    "    \n",
    "The authors put their paper under a Creative Commons license, and [make their data available here](http://discovery.ucl.ac.uk/1526521).  I thank them for their commitment to open, reproducible science, which made this case study possible.\n",
    "\n",
    "The goal of their paper is to model the thermal behavior of a wall as a step toward understanding the \"performance gap between the expected energy use of buildings and their measured energy use\".  The wall they study is identified as the exterior wall of an office building in central London, [not unlike this one](https://www.google.com/maps/@51.5269375,-0.1303666,3a,75y,90h,88.17t/data=!3m6!1e1!3m4!1sAoAXzN0mbGF9acaVEgUdDA!2e0!7i13312!8i6656)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following figure shows the scenario and their model:\n",
    "\n",
    "![Figure 2](https://ars.els-cdn.com/content/image/1-s2.0-S0378778816313056-gr2.jpg)\n",
    "\n",
    "On the interior and exterior surfaces of the wall, they measure temperature and heat flux over a period of three days.  They model the wall using two thermal masses connected to the surfaces, and to each other, by thermal resistors.\n",
    "\n",
    "The primary methodology of the paper is a Bayesian method for inferring the parameters of the system (two thermal masses and three thermal resistances).\n",
    "\n",
    "The primary result is a comparison of two models: the one shown here with two thermal masses, and a simpler model with only one thermal mass.  They find that the two-mass model is able to reproduce the measured fluxes substantially better."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tempting as it is, I will not replicate their method for estimating the parameters.  Rather, I will implement their model and run it with their estimated parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cells download and read the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "download('https://raw.githubusercontent.com/AllenDowney/' +\n",
    "         'ModSim/main/data/DataOWall.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv('DataOWall.csv', \n",
    "                   parse_dates=[0], index_col=0, \n",
    "                   header=0, skiprows=[1,2])\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The index contains Pandas `Timestamp` objects, which is good for dealing with real-world dates and times, but not as good for running the simulations, so I'm going to convert to seconds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "timestamp_0 = data.index[0]\n",
    "timestamp_0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Subtracting the first `Timestamp` yields `Timedelta` objects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_deltas = data.index - timestamp_0\n",
    "time_deltas.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can convert to seconds and replace the index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.index = time_deltas.days * 86400 + time_deltas.seconds\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The timesteps are all 5 minutes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.all(np.diff(data.index) == 300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the measured fluxes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.Q_in.plot(color='C2')\n",
    "data.Q_out.plot(color='C0')\n",
    "decorate(xlabel='Time (s)',\n",
    "         ylabel='Heat flux (W/$m^2$)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the measured temperatures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.T_int.plot(color='C2')\n",
    "data.T_ext.plot(color='C0')\n",
    "decorate(xlabel='Time (s)',\n",
    "         ylabel='Temperature (degC)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Making the System object\n",
    "\n",
    "`params` is a sequence with the [estimated parameters from the paper](https://www.sciencedirect.com/science/article/pii/S0378778816313056#tbl0005)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "R1 = 0.076   # m**2 * K / W,\n",
    "R2 = 0.272   # m**2 * K / W,\n",
    "R3 = 0.078   # m**2 * K / W,\n",
    "C1 = 212900  # J / m**2 / K,\n",
    "C2 = 113100  # J / m**2 / K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = R1, R2, R3, C1, C2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll pass `params` to `make_system`, which computes `init`, packs the parameters into `Series` objects, and computes the interpolation functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_system(params, data):\n",
    "    \"\"\"Makes a System object for the given conditions.\n",
    "    \n",
    "    params: Params object\n",
    "    \n",
    "    returns: System object\n",
    "    \"\"\"\n",
    "    R1, R2, R3, C1, C2 = params\n",
    "    \n",
    "    init = State(T_C1 = 16.11, T_C2 = 15.27)\n",
    "    \n",
    "    t_end = data.index[-1]\n",
    "    \n",
    "    return System(init=init,\n",
    "                  R=(R1, R2, R3),\n",
    "                  C=(C1, C2),\n",
    "                  T_int_func=interpolate(data.T_int),\n",
    "                  T_ext_func=interpolate(data.T_ext),\n",
    "                  t_end=t_end)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make a `System` object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "system = make_system(params, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test the interpolation function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "system.T_ext_func(0), system.T_ext_func(150), system.T_ext_func(300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementing the model\n",
    "\n",
    "Next we need a slope function that takes instantaneous values of the two internal temperatures and computes their time rates of change.\n",
    "\n",
    "The slope function gets called two ways.\n",
    "\n",
    "* When we call it directly, `state` is a `State` object and the values it contains have units.\n",
    "\n",
    "* When `run_solve_ivp` calls it, `state` is an array and the values it contains don't have units.\n",
    "\n",
    "In the second case, we have to apply the units before attempting the computation.  `require_units` applies units if necessary:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function computes the fluxes between the four zones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_flux(t, state, system):\n",
    "    \"\"\"Compute the fluxes between the wall's surfaces and the internal masses.\n",
    "    \n",
    "    state: State with T_C1 and T_C2\n",
    "    t: time in seconds\n",
    "    system: System with interpolated measurements and the R Series\n",
    "    \n",
    "    returns: Series of fluxes\n",
    "    \"\"\"    \n",
    "    # unpack the temperatures\n",
    "    T_C1, T_C2 = state\n",
    "        \n",
    "    # compute a series of temperatures from inside out\n",
    "    T_int = system.T_int_func(t)\n",
    "    T_ext = system.T_ext_func(t)\n",
    "    \n",
    "    T = [T_int, T_C1, T_C2, T_ext]\n",
    "    \n",
    "    # compute differences of adjacent temperatures\n",
    "    T_diff = np.diff(T)\n",
    "\n",
    "    # compute fluxes between adjacent compartments\n",
    "    Q = T_diff / system.R\n",
    "    return Q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can test it like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "compute_flux(0, system.init, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a slope function that computes derivatives of `T_C1` and `T_C2`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slope_func(t, state, system):\n",
    "    \"\"\"Compute derivatives of the state.\n",
    "    \n",
    "    state: position, velocity\n",
    "    t: time\n",
    "    system: System object\n",
    "    \n",
    "    returns: derivatives of y and v\n",
    "    \"\"\"\n",
    "    Q = compute_flux(t, state, system)\n",
    "\n",
    "    # compute the net flux in each node\n",
    "    Q_diff = np.diff(Q)\n",
    "        \n",
    "    # compute the rate of change of temperature\n",
    "    dQdt = Q_diff / system.C\n",
    "    return dQdt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test the slope function with the initial conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "slopes = slope_func(0, system.init, system)\n",
    "slopes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's run the simulation, generating estimates for the time steps in the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "results, details = run_solve_ivp(system, slope_func,\n",
    "                                 t_eval=data.index)\n",
    "details.message"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the results look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "results.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_results(results, data):\n",
    "    data.T_int.plot(color='C2')\n",
    "    results.T_C1.plot(color='C3')\n",
    "    results.T_C2.plot(color='C1')\n",
    "    data.T_ext.plot(color='C0')\n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Temperature (degC)')\n",
    "    \n",
    "plot_results(results, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These results are similar to what's in the paper:\n",
    "\n",
    "![Figure 5](https://ars.els-cdn.com/content/image/1-s2.0-S0378778816313056-gr5.jpg).  \n",
    "\n",
    "To get the estimated fluxes, we have to go through the results and basically do the flux calculation again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def recompute_fluxes(results, system):\n",
    "    \"\"\"Compute fluxes between wall surfaces and internal masses.\n",
    "    \n",
    "    results: Timeframe with T_C1 and T_C2\n",
    "    system: System object\n",
    "    \n",
    "    returns: Timeframe with Q_in and Q_out\n",
    "    \"\"\"\n",
    "    Q_frame = TimeFrame(index=results.index, \n",
    "                        columns=['Q_in', 'Q_out'])\n",
    "    \n",
    "    for t, row in results.iterrows():\n",
    "        Q = compute_flux(t, row, system)\n",
    "        Q_frame.loc[t] = (-Q[0], -Q[2])\n",
    "        \n",
    "    return Q_frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q_frame = recompute_fluxes(results, system)\n",
    "Q_frame.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how the estimates compare to the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_Q_in(frame, data):\n",
    "    frame.Q_in.plot(color='gray')\n",
    "    data.Q_in.plot(color='C2')\n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Heat flux (W/$m^2$)')\n",
    "    \n",
    "plot_Q_in(Q_frame, data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_Q_out(frame, data):\n",
    "    frame.Q_out.plot(color='gray')\n",
    "    data.Q_out.plot(color='C0')\n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Heat flux (W/$m^2$)')\n",
    "    \n",
    "plot_Q_out(Q_frame, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These results are also similar to what's in the paper (the bottom row):\n",
    "\n",
    "![Figure 3](https://ars.els-cdn.com/content/image/1-s2.0-S0378778816313056-gr3.jpg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
